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Abstract 

A direct method for calculating the minimal length of "one-dimensional" Josephson 
junctions is proposed, in which the specific distribution of the magnetic flux retains 
its stability. Since the length of the junctions is a variable quantity, the corresponding 
nonlinear spectral problem as a problem with free boundaries is interpreted. 

The obtained results give us warranty to consider as "long" , every Josephson junc- 
tion in which there exists at least one nontrivial stable distribution of the magnetic flux 
for fixed values of all other parameters. 
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I. Posing the Problem 



It is known that the stationary distributions of the magnetic flux ip(x) in "long" 
(one-dimensional) Josephson junctions (JJ) are solutions of the nonlinear boundary 
value problem (BVP) 

- 9xx + jL>0)siny9 + 7 = , x G (-R, R) , (1) 
<p x (±R) = h B , (2) 

where hs is the external magnetic field alongside the axis y on the junction plane (see 
Fig.la). 

We note that the kind of boundary conditions (fj) , either in presence or absence of 
current 7 in the right side of eq. ([]]), are determined by the geometry of the junction. 
Here we consider simple junctions with overlap geometry ,@ in which the current 7 can 
be approximately considered as a constant. The generalization of our results in cases 
of any other geometry, for example in-line geometry, does not require big efforts. 

We suppose that the given continuous function < jo(x) < 1 describes the varia- 
tions of the Josephson current amplitude, caused by the possible local inhomogeneities 
of the dielectric layer thickness. When jo{x) = 1 the junction is homogeneous. Oth- 
erwise when the junction is inhomogeneous the function jn{x) is usually modelled by 
Dirac 5 - function!'! or its continuous approximations, for example hyperbolic func- 
tions,! splines! etc. At the present work as is in,! we use an isosceles trapezium with 
base fj, (see Fig. 16) as more suitable in physical sense model of inhomogeneity. 

Every solution of the equation ([I]) is simultaneously a stationary solution of the 
perturbed Sin-Gordon equation (SGE) 

Vtt + atpt - f xx + j D (x) sin ip + 7 = , (3) 

where a is coefficient of a resistance. When a > the second term in the above equation 
is dissipative and hence arbitrary distribution (p(x,t) of the magnetic flux as result of 
an energy loss can be "attracted" by some steady distribution (f(x). 

In order to study the stability of some concrete solution <p(x) of the BVP (0),(@) we 
consider the following Sturm-Liouville problem (SLP) 

- ^xx + q( x ) ip — Xi/j, xe{-R,R), (4) 

where q(x) = jn(%) cos<£>(x) is a potential, originated by the solution ip(x), and bound- 
ary conditions of Neumann's type 

il> x {±R) = . (5) 

It is well known that the SLP of such kind contains a counting set of different 
eigenvalues 

X-min = Ao < Ai < A2 < ... A n < ... , 

and every one corresponds to a unique eigenfunction if) n (x) , n = 0,1,2, ... , determi- 
nated by the norm condition 



r-R 

<ip,tp>= ip 2 (x)dx = l. (6) 

J-R 

If the minimal eigenvalue X min > , the respective solution ip(x) of BVP (0),© 
is stable with respect to small time-space perturbations. If the minimal eigenvalue 
Xmin < , this solution is unstable. The eigenvalue \ m in = is a bifurcation point, in 
which the stable solutions of eq.([l|) go to unstable ones and vice versa (for details see!). 

Apart from the space coordinate x, the virtual solutions of the nonlinear BVP 
(|I]), (0) depend also on the physical parameters 7 and "technological" ones fi and 
R , i.e., ip = ip(x,p), where we simply substitute p = {hs, 7, /jl, R}- The varying of 
every of those parameters causes a variation of the distribution ip(x,p) and therefore 
subsequent variations of the potential q(x,p), the eigenvalues A(p) and the respective 
eigenfunctions ip(x,p). Thus we can conclude that every solution of BVP ([[]), (Q) has 
an area where it remains stable with regard to the variations of the parameters p. The 
equation 

X m in{P) =0 (7) 

determinates in the parametric space a hypersurface which points appear to be bifurca- 
tion points corresponding to the solution under consideration ([[]). The intersections of 
the bifurcation hypersurface (0), when there are fixed pairs of the parameters p, we call 
bifurcation curves and respective values of the parameters - bifurcation (critical) param- 
eters. The most interesting from the physical viewpoint seem to be bifurcation curves 
"external current - magnetic flux" A m j n (/iB, 7) — 0, when the geometrical parameters \x 
and R are fixed. That curves could be relatively easy obtained experimentally.! 

From the technological point of view, however, it is worth investigating the bifurca- 
tion curves as functions with respect to at least one of geometrical parameters 

Xmin(h B , R) = , Xmin(l, R) = , Or A mirl (/i, R) = . 

Such kind of problems can be connected with the optimization of sizes of devices, 
containing Josephson junctions. 

Formally the numerical modelling of the bifurcation curves as function of some 
concrete parameter po 6 p (at the paperS it is chosen po = hs) can be schemed as 
follows. We find some solution of the BVP (fl|), (0). After that we check his stability 
with respect to small perturbations solving the corresponding SLP (§)-©. If X m in > 
we set a new value po : = po + Ap, where Ap is given increment, and solve the BVP ([[]), 
(0) again using the obtained solution as initial approximation. We repeat this iteration 
while at the "current stage" po the equation X min = is satisfied numerically. Then 
one point at the searched bifurcation curve is calculated. 

At Fig. 2 the typical relationships A mi „(A) (here A = 2R is the junction length) 
corresponding to the so called "main fluxon"[| (see Fig. 3) in inhomogeneous junction 
containing one resistance inhomogeneity placed at the point x = when \i = 0.5 and 
\i — 1 respectively are shown. The points B and B\ appear to be bifurcation points of 

infinite JJ with one <5-shaped microinhomogeneity at the point x = "main" fluxon/antifluxon 
is represented by the exact solution ip(x) — 4 arctanexp(±x)Jj 



corresponding relationships A TO j„(A) = 0. According to the above mentioned reasonings 
these points determine the minimum length of JJ, for which main ffuxon is still stable. 

Obviously such kind of algorithm is quite hard. Therefore it is quite natural to put 
the question how to calculate directly the bifurcation curves. A general approach study 
for the posed problem is proposed at the articles.!' We consider the equations (|J), 
(0), (||) - (Q) as a closed nonlinear system with respect to the functions f(x), i/j(x) and 
one of parameters p (for example, /i^ or ")H), while the other 3 parameters, and A also, 
are given. Then fixing A to be small enough, (for example A = 0.01), every solution of 

d A 

the above system with a priori prescribed accuracy (the derivative 7: ► 00 when po 

opo 

approaches the its critical value) belongs to the small vicinity of the searched bifurcation 
curve. 

The calculation of the critical half-length R of homogeneous or inhomogeneous JJ, 
corresponding to a concrete nontrivial distribution of the magnetic flux, is an important 
practical problem. The shortcoming in this connection ensues from the fact that the 
equations ([]]), (0), ([|) - flf]) are implicit with respect to the quantity R. At the present 
work we propose how to overcome the mentioned imperfection. 



II. Method of Solution 



For given values of parameters A, n, hs and 7 we consider the system ([[]), (0), (U) - 
(H) as a nonlinear eigenvalue problem with respect to the eigenfunctions ip(x),i])(x) and 
to the eigenvalue R. As the parameter R does not occur explicitly , it is convenient to 
use the Landau transformation £ = —^-—.0 Choosing simply f(R) = R we map the 

did 

original interval [—R,R] to the interval [—1, 1]. Taking into account that — = — — 
we obtain that the above system renders to the system: 



-^ + i? 2 [J D (Osin^ + 7 ] = 0, (8) 

tp 6 (±l) - Rh B = , (9) 
+ R 2 [ j D (0 cos ^ - A] -0 = , (10) 
^(±1) = 0, (11) 

<^,i> > -1 = rJ ^ 2 (C)di-i = 0, (12) 

where (7) denotes the function of £. Without fear of confusion the bars will be omitted 
henceforth. Let us consider the left-hand side of the system (H)-(I2) as a functional 
vector F(ip, if), R). Then we have F(<p, ijj, R) = 0. 

The nonlinear system under consideration will be solved using the continuous analog 
of Newton's method.i Following it, we introduce a continuous parameter ("time") 
t G [0,oo) and suppose the quantities ip, ifj and R depend on t, i.e., ip(t,£), i>(t,£,) and 
R(t). The following abstract differential equation is used 

dF 

— + F = F^ + F^ + F' R R + F = 0. 



Here 



d 



F' ip u= -F(<p + eu) 



, K v 



e=0 



d 

de 



F(ip + ev) 



e=0 



,F> R p=jF{R+ep) 



t=0 



are Frechet's derivatives of F, and 



u — <p, v = if), p = R 



(13) 



are the "time" derivatives of the functions <p, tp and R, respectively. Simplifying above 
system we obtain 



-m k + R 2 j D (0 cos up u + {2R \j D (g) sin ip + 7] + R 2 j D ,z(0 sin ip} p 

-^ + R 2 [jd(0^^ + i} = 0, 

li C (±l) - + ^(±1) - i2fc B = , 

-VX + R 2 [q(£)-\]v+{2R [q(0-M^ + R 2 joAO^cosip} p 
-R 2 j D (0 iisincpu- + R 2 [g(0 - A] ^ = , 
^(±1)+^(±1) = 0, 

2R < ip,v > + p < i),if) > + R < i/j,i> > -1 = , 

This system can be solved using the following decomposition: 

u = Ui + p u 2 , v = Vi + p v 2 , 

where Wi(0> u 2(0> ^(0 are new unknown functions of £. That assumption yields 

four linear two-point boundary problems with respect to the new introduced functions: 



(14) 
(15) 

(16) 
(17) 
(18) 



-u m + R 2 q{0 Ul = ^(0 - R 2 []d{0 sin^(0 + 7] 
Ul ^±l) = Rh B -(p ( (±l), 



(19) 



-u^ + R 2 q(0 u 2 = -2R [j D (0 sin <p(£) + 7] - R 2 (£) £ sin <p(£) 

tt2f(±l) = , 



(20) 



+ ^ [9(0 - A] «i = - ^ [7d(0 008^(0 - A] m + 

R 2 j D ,dom^^(OMO 

v u (±l) = -^(±1) , 



(21) 



-v 2 ^ + R 2 [q(0 - A] v 2 = R 2 JD (0 m sin^(e) u 2 (0 - 2R [<?(£) - A] ' 
R 2 3dA0^(0 cosp(0 

(22) 

v 2f (±l) = . 
At last the norm condition renders to 

1 - R < ip,ip > -2R < ip,v x > 

P— • (23) 

' < ip, ip > +2R < ip, v 2 > 



For numerical computation the time derivatives u, v and p at (|T3| ) are discretizated 
by Euler's method (see details ini). For given iteration k, the approximation at the 
next stage k + 1 is obtained as follows 

R k+i = R k + T k p k ^ ^fe+i = ^k + T k u k ^ ^k+i = ^k + r k v k_ 

Here r fc e (0,1] is a parameter ("time" step) which can be chosen satisfying the 
condition the residual to be minimal (seel). 

Let us note that every one of BVP (19) -(22) can be presented simply in the form 

-y& + p(€)y = r(£), 

(24) 

%(±i) = y±, 

where p(£) , r (£) are given functions, and y± are given constants. Further we define 

an uniform set at the interval [—1,1] namely £„• = — 1 + jh , j = 0, . . . , N, where N 

2 

is number of knots, h — — is the step of set. Let £>(£) is a cubic spline interpolating 

the function ?/(£) over the £-mesh. We assume that M(£) =: S'^(^). Then taking into 
account the continuity condition concerning the first moments of the spline S(£) and 
BVP (PJ), we obtain a three-diagonal algebraic system (se 

The iteration process starts from the initial conditions ipj = (fj, ip® = ipj, j = 1, . . . N 
and R° = R k , where k denotes the number of the iteration. We solve consequently the 
two-point BVP (19) -(22). Thus the values for the grid functions u\, ,u\^ , v\, are 
obtained. Then using (^) we calculate the increment p. Hereupon we calculate the 
increments u and v and obtain the predictions for the junction length R k+l and the grid 
functions , at the new stage fc + 1. The criterion for terminating the iteration 
is 



max (S(p, Sip, 5R) < e , 



where S(p, dip and 5R are the corresponding residuals. We use the norm estimation 
e ~ 10- 8 -T- 10- 12 . 



III. Results and Discussion 



The numerical correctness of the used scheme is verified through appropriate nu- 
merical experiments. We used different meshes with sizes N = 256, 512 and 1024. 
The relative differences for the magnetic flux (p(x), first eigenfunction ip(x) do not ex- 
ceed 0.004% and 0.03%, respectively. We estimated the order of approximation of the 
obtained solution using the Runge method. The calculations carried out on meshes 

with spacings h = — — , — —, — — are presented on the table. It is easy to show the 

±jli(j .ZOO D.L.Z 

approximate relationship 

Zh - Zh 

2 - « 2 2 



Zh — Zh 



holds at the knots 0; y; N. Here z = {cp, ip, R} . On the ground of that comparison we 
conclude the second-order approximation for the functions <p(x), ift(x) and the junction 
half-length R is satisfied. 



N 


h 


tp = tp(-R) 


IfN = <p(0) 


<p N = ip(R) 


% = = 


tpN = V>(0) 




256 


i 

128 


1.1770473 


3.1415927 


5.1061380 


0.4203551 


0.5178962 




512 


1 

2.J6 


1.1770277 


3.1415927 


5.1061576 


0.4202415 


0.5177725 




1024 


5l2 


1.1770244 


3.1415927 


5.1061609 


0.4202198 


0.5177485 





All results, stated below, are related to the solutions of kind "main fluxon" (see Fig. 

3). 

The upper curve on the Fig.4 presents the initial distribution of the magnetic field 
(p x (x) alongside the junction. This distribution is chosen to be a solution of BVP ([[]), 
(0) for R — 5, J%b — and 7 = 0. At the same figure the lower curve presents the 
calculated distribution of the magnetic field corresponding to the minimal eigenvalue 
Ami™ = 0.01. Hence this is the searched bifurcation distribution and in the framework 
of this model the value of the spectral parameter R m i n ~ 2.11 is the minimal half-length 
of the junction providing a stable main fluxon. In this sense, junctions, whose length 
lies below the critical R m in, it is necessary to be considered as "short" for distributions 
of the kind "main fluxon". As it was shown,! for such short length, there is a unique 
stable distribution of the magnetic flux in the junction - Meissner's solution^ 

On Fig. 5 the distributions of magnetic field <p x (x) alongside the junction in absence 
of external current (7 = 0) depending on its boundary values hs = and hs — 1 
respectively, are represented. It is seen that two solutions could be approximated by 
means of two-degree polynomials. Furthermore they seem to be geometrically similar. 
Namely the curve <p x (x) corresponding to magnetic field hs > may be considered as 
obtained from the other curve doing a homothety alongside axis x with respect to the 
pole x = and a translation alongside the axis <p x (x). 

On Fig. 6 the obtained relationship A m j n (/iB,7) = is drawn, depending on the 
length of the trapezium base fi. It is well noticeable the stabilizing influence of the 

2 This is the trivial solution <p(x) = (stable) and y?(x) = n (unstable) for hs = 
and 7 = 0. 



boundary magnetic field upon the critical length of the junction, i.e., increasing the 
magnitude of above mentioned quantity one decreases the critical length of junction 
2R min . Thus the calculated critical half-length of the Josephson junction corresponding 
to the so called main soliton (Jib = 0, 7 = 0) is R m in ~ 2.11 while the same quantity 
when h B = 1 and 7 = is R m in ~ 1-35. 

IV. Concluding Remarks 

A direct iterative method for obtaining the minimal half-length R m i n , corresponding 
to fixed distribution of the magnetic flux in a long Josephson junction is developed. The 
mere existence of the minimal length is a warranty enough for us to name "long" , every 
JJ in which there exists at least one nontrivial stable distribution of the magnetic flux 
for given values of other parameters. 

An appropriate linearization based on the continuous analog of Newton's method 
renders the original nonlinear spectral problem to four two-point linear BVP. A spline- 
difference scheme in second order of approximation for solving of these BVP is used. 

This method may be applied for solving more general nonlinear problems. 
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VI. FIGURE CAPTIONS 

Figure 1: a) Geometrical sketch of an inhomogeneous J J; b) Geometrical model of 
the amplitude of Josephson current. 

Figure 2: The minimal eigenvalue X m i n as function of the junction length A = 2R. 



Figure 3: The magnetic flux (p(x), magnetic field (fi x (x) and first eigenfunction tp(x) 
as functions of the junction length A = 2R. 

Figure 4: The magnetic field <p x (x) alongside the junction as function of its length 
2R: "V" - initial distribution; "o" - final (bifurcation) distribution, corresponding to 
the minimal length. 

Figure 5: The magnetic field <f x (x) alongside the junction as function of its length 
2R for length of the trapezium base n — 1, bias current 7 = and different values of 
the boundary magnetic field Kb'- " o " - Kb — ; " o " - hs — 1 ■ 

Figure 6: The minimal (bifurcation) eigenvalue \ m in as a function of the boundary 
magnetic field Jib for different lengths of the trapezium base [x: "V" - \x = 0.5 ; "A" - 
H = l. 
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